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deformation; (2) gravity variation; (3) deflection of the vertical; and (4) strain tensor which 
is induced by these deformations. The closed form sum of all those harmonic series which 
are required for such numerical evaluations are obtained and displayed in tabular form. 

We finally examine the liquid core/solid mantle interaction by studying the asymptotic 
expansions of the equations of motion as the frequency of the forcing function tends to zero. 
We characterize such interaction as being generated by a boundary layer. Such a layer 
justifies the discontinuity that certain variables must undergo across this interface in order 
to numerically integrate the Navier-Stokes equations with boundary conditions. 
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LOAD NUMBERS, SOLID EARTH TIDES, 
AND LIQUID CORE DYNAMICS 


INTRODUCTION 

In NRL Report 8410 [1] we expressed the deformations of an elastic Earth due to static surface 
loads as infinite series of spherical harmonics. The coefficients of these harmonics were called the load 
numbers and were obtained by integrating the Navier-Stokes equation throughout the interior of a lay¬ 
ered Earth with boundary conditions that had to be satisfied at both the center of the configuration and 
the loaded upper surface. To numerically evaluate these infinite series, one must: (a) be able to ascer¬ 
tain the load numbers of arbitrarily high order without having to solve each time a boundary problem 
(this is tantamount to establishing an asymptotic expansion for the sequence of load numbers) and (b) 
develop practical and efficient ways for numerically summing those series of spherical harmonics once 
the load numbers are known. 

We discuss in this report a question that was left unanswered in our previous publication, that is: 
why the stress function was assumed to be discontinuous at the core-mantle interface when numerically 
evaluating the load numbers as a boundary value problem. 

This report presents in a systematic and rigorous way procedures suitable for the numerical 
evaluation of the^ Earth’s sph eroidal deformations, which most commonly are called the Earth’s tides. 
For this purpose, we: (a) study the Boussinesq theory for the elastic deformations of a flat plate and 
show how by applying the Boussinesq solution to the tangent plane at a loading point of the spheroidal 
Earth we can establish the asymptotic behavior of the load numbers [2]; (b) elaborate on certain 
numerical procedures which seem to be best suited for the summation of series once the asymptotic 
behavior of their coefficients has been ascertained; (c) provide in tabular form the closed form expres¬ 
sions of those infinite series of spherical harmonics which are essential for the numerical evaluation of 
the Earth’s tides; and (d) examine the dynamics of a liquid core for a nonrotating Earth to mathemati¬ 
cally prove the existence of a boundary layer at the top of the liquid core. This layer provides a 
justification for the discontinuous behavior of some of the integration variables at the liquid core/solid 
mantle interface, and furnishes the proper number of free parameters for satisfying the boundary condi¬ 
tions at the loaded surface. 

/K 

A more detailed discussion of these ideas can be found in a forthcoming publication by Lanzano 
[31. In reaching our main results, we want to acknowledge the fact that we have greatly benefited from 
two older publications by Farrell [4] and Pekeris and Accad [51. 

ASYMPTOTIC VALUES OF THE LOAD NUMBERS 

The strategy we follow in order to ascertain the asymptotic values of the load numbers can be 
summarized in the following steps: 

(I) One should compare the infinite series expansions which represent the spheroidal defor- 
. mations and which contain the load numbers as the coefficients of the various harmonics 
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with appropriate closed form solutions for the same type of deformations which are valid 
within a limited region of the spheroid. This is the case in the neighborhood of a loading 
point where the Boussinesq theory of the flat plate applies and gives rise to a closed form 
solution. 

(2) One should refer both solutions in the region where they overlap to a common reference 
frame. It is convenient to make use of the cylindrical coordinate system with the loading 
line as its symmetry axis and the tangent plane to the spheroid at the loading point as the 
auxiliary plane. 

(3) Because the limit of the Legendre polynomials P n for large values of n is the Bessel func¬ 
tion Jo (of order zero) [6], and similarly the same limit for the derivative dPjd9 is pro¬ 
portional to J { , one should solve the Boussinesq problem by means of the Hankel 
transforms of order zero and order one because these functional transforms use those 
Bessel functions as their kernels. 

(4) As a consequence, it follows that the asymptotic values of the three load numbers must 
be proportional to the Hankel transforms of the two components of the displacement and 
of the perturbed potential in the Boussinesq problem. 

We shall next elaborate statements 1 through 4. To begin with, it is clear that for small angular 
separations from the loading line, the displacements and the perturbed potential of a layered spheroid 
must agree with the corresponding quantities referring to the tangent plane at the loading point. In the 
neighborhood of a loading point, and for the values of density and elastic parameters commonly 
accepted in the most recent Earth models, the elastic forces dominate the gravitational forces. We can, 
therefore, not only assume that the density and Lamd parameters remain constant in the neighborhood 
of a loading point, but also neglect the coupling between the gravitational and elastic forces; this situa¬ 
tion enables us to study and separately solve one elastic problem and one gravitational problem. 

In the absence of the gravity field g 0 of the reference state and of the perturbed gravity gj due to 
deformation, assuming also: (1) elastic equilibrium, and (2) constancy of the Lam6 parameters, the 
Navier-Stokes equation, which was developed in our previous work, Eqs. (14) and (15) of NRL Report 
8410 [1], will simply reduce to 

(X + 2/4.)V(V • V) - mV x V x V, (1) 

where u is here the displacement vector. The above equation should be solved for the material half¬ 
space simulating the tangent plane, under the assumption that there is only one normal load concentrat¬ 
ed at a given location in the uppermost surface. By the same token, the perturbations in the gravita¬ 
tional potential can be represented by the Poisson equation 

V 2 <*> - -4irGp„(V ■ 17). (2) 

with constant density p 0 , along the material half-space, and by the Laplace equation 

V 2 * - 0, 

in empty space. 

Consider a spheroidal Earth bounded by an equipotential surface_r — a with a vertical load applied 
at one of its points Q. We introduce the spherical reference frame (OQ, r, 9) with origin Oat the cen¬ 
troid of the Earth and polar axis along the line OQ, we shall measure the^olatitude 9 from this axis and 
assume that the configuration be symmetrical with respect to the line OQ. This line of application of 
the load can be taken as the r-axis of a cylindrical reference with origin at the loading point Q, and the 
tangent plane to the spheroid at Q can be taken as the auxiliary plane upon which the other two vari¬ 
ables (rand *) will be measured. This cylindrical system shall be denoted by (Q, r, <li). It is clear that 
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if we take a point T in the neighborhood of Q, the arc (QT) along the equipotential, which is 
represented by a9 , must be the radial distance in the cylindrical reference system; also that the unit 
vectors e z and e r of the cylindrical frame must be the limiting positions of the unit vectors e, and e», 
respectively, of the associated spherical frame. 

We shall solve both Eqs. (1) and (2) where the displacement vector V will have the following 
representation in cylindrical coordinates: 

V - u(, 2 ,r)e z + v(z,r)e r . (3) 

The spheroidal deformations of a layered spheroid along an equipotential surface r — a due to an 
applied vertical load of mass mo acting at a point Q, and when expressed in spherical coordinates in 
terms of the load numbers h\ /', can be written as follows: 

Via,9) = T [e, h'„{a)P„(cos 9) + e e l‘„{a)dP„{cos 9)/d9], (4) 

m n-0 

Here m is the mass of the spheroid within the equipotential surface r — a. 

Similarly, the potential due to the deformation of the spheroid can be expressed in terms of the 
third load number k\ also in spherical coordinates, as follows: 

<t>(a,9) -- y 0 (a) T k„(a)P n ( cos 9). (5) 

m n-0 

Here yo(a) is the gravitational attraction since we are neglecting the rotational effects. 


It is now appropriate to recall that the limit of the Legendre polynomials for large values of the 
parameter n can be established in terms of the Bessel functions. More specifically, see Ref. 6, p. 155 


lim P n [cos (9/n)] - J 0 (9). (6) 

IT—oo 


By differentiation, and if use is made of elementary properties of Bessel functions, 

-J,(9). 


1 dP„ [cos (0/n)] 
-I™ n dWn ) 


(7) 


Because of these properties it is convenient to solve Eqs. (1) and (2) by means of the Hankel 
transforms. 


The Hankel transform of order n of a function f(z, r) of two variables with respect to one of its 
variables is represented by the corresponding capital letter and is defined as 

F(z, *) - f Q fiz, r)J„{(r)rdr, (8) 

where J„ is the Bessel function of order n. We can verify that the inversion of Eq. (8) will provide 

/fc r)-f~F(z. ()JMr)(d(. (9) 

Note that the product (r of these two variables must be a nondimensional quantity. From Eqs. 
(8) and (9), it is clear that the dimensions of a set of Hankel transforms (/;/") will differ by the square 
of the dimension of the f or r-variable. 

We now revert to Eq. (1) and write the components of the stress as 

r a - A(e„ + < M + f B ) + 2pf B - (A + + ~ 

He 
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Here u , v are the components of the displacement vector 17given by Eq. (3); the notation we use is the 
one by Sokolnikoff [7]. 

Consider the system formed by Eqs. (1) and (10); and apply the Hankel transform of order zero 
to the u and z a variables, the H transform of order one to the v and r a variables. Denote the H 
transforms of these four variables by the corresponding capital letters. Then assume that both u and v 
and their partial derivatives vanish at infinity with Mr. Consider the system of equations in the 
material region of the linear space, which we assume to be the z < 0 half space. 

We reach the linear differential system 

T- AS, 

where the primes denote derivatives with respect to the ^variable and where we have set 

3? s Column (U, V, T a , T a ). 

The matrix A is 

0 —Xf/<r 1/cr 0 

* o oi/m 
0 0 0 -£ 

0 A/j.r)^/cr Xf/o- 0 

with constant values for <r - X + 2 m and tj - X + ft. 

The general solution of Eq. (11) is obtainable by elementary procedures in terms of the eigen¬ 
values p of the A matrix, i.e., in terms of the roots of the characteristic equation 

Det 04 - pi) - 0, (14) 

where I here stands for the (4 x 4) identity matrix. The four roots of Eq. (14) are 0i - Pi - f and 
03 - 04 --*. 




( 11 ) 

( 12 ) \ 


The eigenvectors or solutions to Eq. (11) can be expressed according to the infinite relationship 


3? - exp (±fz) exp [(A =F (I)z]y - exp (±fz)l? + z(A =F f I)y + y(/4 =F f/) 2 ? + . 

..]. (15) 


One can obtain two independent solutions. 


"V 

exp (±{ z)7 x , 

by solving the linear system. 

(16) 


(A T £/)?, - 0. 

(17) 


Two additional independent solutions are 


l 

exp (±{z)17 2 + zC 4 Tf/)y 2 ], 

(18) 

i 

obtainable by choosing solutions of 


i 

i 

U T (I) 1 ?! - 0, 

(19) 

) 


which are not common to Eq. (17). In both cases, the infinite expansion appearing in Eq. (IS) will ter- 
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By following the preceding procedure, we reach the set of four linearly independent solutions: 


1 

=fl 

±2 fi£ 

-2#if 


exp (±fz); 


( 20 ) 


Tfz +n/ V 

fz ± o7t) 
"2 H( 2 z 
±2n$ 2 z + 2 fi( 


exp (±fz). 


Let us now take into account the boundary conditions. The vanishing of the four original variables at 
z — —oo, means that the transformed variables U, V, T a , T n must vanish at minus infinity. If we 
choose henceforth f > 0, then the solutions containing exp (-fz) cannot be taken into consideration. 
To combine the remaining two solutions, we suppose: (a) the transversal component of the stress van¬ 
ishes at z — 0: r n (0, r) — 0; and (2) the unit normal load is applied on a spherical cap of radius r — 
a, i.e., 

r a (0, r) »> — 1/w a 2 for r < a, 
r a (0, r) - 0 for r > a. 

It follows then that 


and 


TJ0.&-0. 


T b (0,() - l — C J Q {ri)r dr 

it a * u 


1 2 J\(£a) 

2l r fa 


whose limit for a shrinking a is — — l/2ir. Using the previous relations, we finally have 


U(z.O 


1 



exp (fz); 


K(z,f) 


1 


fz + exp (fz). 


( 21 ) 


We next consider the gravitational problem by solving the Poisson equation within the material tangent 
plane (z < 0) and the Laplace equation in empty space (z > 0). 

We introduce the Hankel transform of order zero for the perturbed potential 4> and denote it by 
d>. The Poisson equation then becomes 



<Mz,f) 


2<jpo 

V 


exp (fz). 


( 22 ) 
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after use has been made of Eq. (21) to express the components U and V of the displacement. The 
characteristic exponents of Eq. (22) are ±£. The solution of the Laplace equation valid for empty 
space (z > 0) must vanish at z = +°° and must therefore be in the form 

<J>,(z,f) — A exp (— |z). (23) 


The solution of Eq. (22), valid within the elastic medium, can be seen to be: 




B + 


Gpq | 

vi 1 


exp (fz). 


(24) 


However, must be continuous at the interface z — 0; this leads to A — B. By rewriting the Poisson 
equation as the divergence of the vector 


V0 + 4w GpoV, 

we realize that the normal component of such vector must be continuous at z - 0. In terms of the H 
transforms, we reach the condition: 


ld4> e /dz lj_o” [d<J\/3z] 2 _o + AnGpoUjiO.g): 


(25) 


this is so because of the vanishing of U t (z,£). Use of Eq. (25) is instrumental in determining the 
value of A. We get 


<Mz.f) 


Gp(\ 

- \ exp (-fz) for z > 0, and 


*/Uf) 


Cpo 

2pi 2 


1 + 



exp (f z) for z < 0. 


(26) 


We are now in the position of comparing two solutions which are both valid in the neighborhood of a 
loading point and of determining the limit of the three load numbers as the order of the harmonic 
tends to infinity. 


First, by expressing the solutions given by Eqs. (4) and (5) in cylindrical coordinates and by mak¬ 
ing use of Eqs. (6) and (7), the displacement vector and the perturbed potential for large values of n, 
may be written: 


—- le z h'(a)Jo(n9) - e,nl'„(a)J\(n9))\ - - y a (a)k n (a)J 0 (nO). (27) 

m m 

We next represent the Boussinesq solution for z = 0 and for small values of the radial distance r — a9 
as integrals of the H transforms according to the following relationships: 


S 1(0,a9) - e : f“ i/jo.-j] /<,(«•)-jr dn + e, J" pjo.-jJ./, (»*)—- dn\ 
*(O,a0) - J 0 ” $ 0,^-j J o (,n0)-^ dn. 


(28) 


Note that the integration variable f has been relabeled f — n/o, a legitimate replacement since f has 
the dimension of an inverse length. 
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Each integral appearing in Eq. (28) can be visualized as an infinite summation; their terms per¬ 
taining to large values of n must coincide with the corresponding terms appearing in Eq. (27). Since 
the Boussinesq solution was written in terms of a unit force whereas the applied load here is m a y 0 (a), 
we reach the following relationships: 

hn(a) “ P" 

—- nl'„(a) - ~ -Tf I'jo,— I m 0 y 0 (a)\ (29) 

—~ y 0 (a)k'„(a) - — y 4>|o,— |m 0 yo(a). 
m a 1 \ aj 

By using Eqs. (21) for U and V, and Eqs. (26) for <!>,, we finally arrive at the asymptotic representation, 


A* 

/* 

myo 

1 

A* 

Aira 1 ^ 

~ 3p 0 T|/2jx,p 


Here a — A + 2/a, 17 — X + /x, p is the mean density, m is the total mass of the spheroid, and A*, /*, k * 
are the limits of h„, nl'„, nk'„ as n approaches infinity. All quantities appearing in Eq. (30) must be 
evaluated at the outermost surface, r — a. 

For the 1066A Earth model, these limits can be calculated and are 

A* ** — 11.35767, 

/* - 3.43096, (31) 

A* « -4.70473. 

NUMERICAL EVALUATION OF THE GEOPHYSICAL PARAMETERS 

Once the asymptotic behavior of the three load numbers has been ascertained, one can numeri¬ 
cally evaluate the elements of geophysical interest; these are: (a) the elastic deformations of the sur¬ 
face; (b) the variation of the gravity field in direction and intensity; and (c) the strain tensor induced 
by the deformations. All these quantities are represented by infinite series of Legendre polynomials, 
the angular distance 9 being measured from the loading line. 

Because we know the asymptotic limits A*, A*, /* of h‘„, nk „, we can take advantage of the 
Kummer method for the summation of an infinite series of functions 18,9). This method consists of 
subtracting from the given series an appropriately chosen series of known sum whose elements are 
asymptotically proportional to the series in question. The result can then be expressed as the sum of 
two infinite series of Legendre polynomials: (a) the first one is multiplied by the known asymptotic 
value of the load number and its sum can be represented in closed form; (b) the second infinite series 
has coefficients asymptotically tending to zero, so that it can ultimately be evaluated as a finite sum; the 
number N of its terms must be so chosen that for values of the subscript n larger than N, the difference 
between the corresponding load number and its asymptotic value can be considered negligible according 
to the degree of approximation we plan to achieve. 


Let us further examine these ideas. Consider the components u, v of the displacement vector 
V (a,6). For the normal component, we get 


uUO)i kHa)P„ B (a* i Pn + £ <*»’ - h ' )P * ■ 


(32) 


the value of N depending on the accuracy to be attained. 
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In a similar manner, the tangential component v of the displacement can be evaluated according 
to the scheme: 


v(a,0) 


am 0 


/* T — 

n 


I 

n d9 


I «' 


/•) 1 ^ 


(33) 


This is so because /* is the limit of «/„' as n goes to infinity. Both infinite series, P„ and 
(1/n) ( dPJdQ ), that appear in Eqs. (32) and (33) can be summed in closed form; their expressions are 
given in Table 1. 


Table 1 — Closed Form Sum of Certain Harmonic 
Series of Geophysical Interest 


(1) V P „(cos 0) - y cosec ~ 
n-o 2 2 


(2) T nP„ (cos 0) = - -j cosec I 

-Tn 4 | 2 

,2 


,31 


cosec 


(4) £ 

n— 1 
oo 

(5) £ 

f»-I 

oo 

( 6 ) £ 

1 


I 

n dO 

1 ^ 


cos 9 + sin |y 


cosec0 


n d9 2 

— P„( cos 0) 
n 


1 - sin 3 ^1| cosec 2 0 


-In 


sin 

9_ 

... 0 

1 + sin — 


2 

2 


Next, consider the variation in the gravity field. The potential at a point Q(a + u, 9) of the 
deformed surface, located at a height u from the point P(a,9) on the equipotential r — a, can be writ¬ 
ten: 

Y(Q) = F 0 (P) + uVo(P) + K*(/>) + V'(P). 

Here, primes denote derivatives with respect to the normal distance, V 0 is the potential of the original 
configuration, F* is the potential due to the redistribution of mass, and V, is the potential of the 
applied load. 

A gravimeter located at Q measures the following acceleration: 

V(Q) - Vq ( p) + uVq(P) + ( V’)'{P ) + y;(P). (34) 

We now expand both F* and V, into spherical harmonics and note that Vq (P) — -yo (a), where 

y 0 (a) - J 0 p 0 (r)r 2 dr. (35) 

By differentiation, we easily establish that 

yo (a) " — — y<>(a) + 47 rGpo(fl). (36) 

a 
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If we neglect the contribution due to the density at P , Eq. (34) can be rewritten as 

-y(0) - -y 0 (P) + ~ uy 0 (P) - — V' X (P) + - K,(/>), 
a a a 

which leads to 

y'iP.Q) - y 0 (F) “ y<<?) - — y 0 (/*> £ [2*;(a) - (a + 1)*>) + «]!»,. (37) 

Finally, the tilt or deviation of the vertical direction can be expressed as the angle between the normal 
to equipotential surface, V 0 + V‘ — constant and the geometrical normal to the deformed surface. For 
simplicity, we shall limit our considerations to the component of the tilt upon the meridian of the 
spheroid which we denote by /* (a, 0). 


The normal line to the equipotential surface is represented by 

(1 + O(dPjdO), 

whereas the normal to the deformed surface depends on 

(1 /a)(dr/dO), 

which is proportional to 

hi(dPJde). 

We can then write 


t^a.O) - — f (1 + *„’(«) - A») (38) 

m do 

The numerical evaluation of Eqs. (37) and (38) requires the sum of three additional infinite series. 
nP „; (1 /n)P„\ dPjdO. Their sums are also available from Table 1. 


We now come to the evaluation of the components of the strain tensor. We use the fundamental 
Navier-Stokes equation which was expressed as Eq. (39) in Ref. 1 to write: 


«„(a, 0 ) » u’fa.O) — £ U„(a)P n 
0 


I 

0 


2A (a) ,r / \ i / , \(u) 

UAa) + n(n + 1) - 


r (a) 


a<r(a) 


W + 


E„(a) 


r (o) 


P. (cos 0). 


As usual, primes denote derivatives with respect to the radial distance and cr — A + 2/u.. U, V and E 
denote the expansions into spherical harmonics of the components «, v, of the displacement vector and 
of the component of the stress tensor. We know that E„(a) is a delta-function with its peak at 0 - 
0, so that it will give a zero contribution anywhere else. If we take these facts into account and 
represent the v component of the displacement vector in terms of the load numbers, we can rewrite the 
previous equation as 


e " (a ’ 9) - ~ u(a ’ 0) + ? i nin + l) '»' (fl)/> « (cos 9) - 


We can next verify that 

1 

- 


u(a,0) + — v(a,l 0 ) 


1 , , Wn “ ,, , d*P. 

- S % <■“ ST- 


(39) 


(40) 


Also 




— lu(a, 0) + v(a, 0 ) cot 0], 
a 


9 


( 41 ) 






PAOLO LANZANO 

The other components ot the strain tensor can be shown to vanish when evaluated at the surface r - a. 
The numerical computation of these infinite series which represent the components of the strain tensor 
necessitates the knowledge of the sum of the series (1 / n){d 2 Pj d0 2 )\ t’ : s is also given in Table 1. 

Table 1 provides the closed form sum of all those series which are required for the numerical 
evaluation of the geophysical parameters. A detailed derivation of the sum of these series appears in 
Ref. 3; for some partial results see Ref. 10. 

From a computational point of view, few remarks are in order before bringing this section to a 

close. 


(1) The second infinite series that was introduced through the Kummer transformation turns out to 
be slowly convergent; various numerical algorithms must be used in order to accelerate its conver¬ 
gence; in this regard, the Euler transformation was found to be a useful tool [8,9]. 

(2) For large values of the order n of the harmonic and for small values of the angle 0 measured from 
the loading line, one can achieve computational economy by approximating P„(cos 0) and dPJdB 
by means of J 0 and J\, as has been described in the previous section. 

(3) The functions appearing in Table 1 are not finite at 8 — 0 and should not be used in the immedi¬ 
ate neighborhood of the loading position. We have remedied to this situation by employing the 
Boussinesq approximation extended to a suitable neighborhood of the loading point. We resume 
therefore our considerations on the Boussinesq solution which had reached the representation 
given by Eq. (21) and try to express it in terms of cylindrical coordinates at the loading point. In 
other words, we shall invert the Hankel transformation. For this purpose, we use the following 
integral appearing in Ref. 6: 

J 0 exp(-<jz) J y (bz)dz «• (C/bYia 2 + b 2 )~' /2 , 
which is valid for v > 0, a 0, b ^ 0 and where C — ( a 2 +b 2 ) i/2 - a > 0. 


Differentiation of the above expression with respect to the parameter a yields 

J" fl exp (— az)J y (bz)z dz — [pb(C/b) v ~ l + a(. 1 — p)(C/b)*](.a 1 + b 2 )~ 2/2 . 

The preceding two integral relations for v — 0, I are all that is required to obtain the Boussinesq 
approximation in terms of cylindrical coordinates (z, r): 


u(z, r) 


m 0 y 0 (a) 

4ir/iR 



(42) 


v(z, r) 


m oyp (a) . j_ t )r 2 z 

4 irrjr R M R 3 1 


Here we have set 


R 2 — r 2 + z 2 ; rj — X + n and <r — X + 2 fi. 


LIQUID CORE/SOLID MANTLE DYNAMICS 

To justify certain assumptions that were made in carrying out the numerical integration of the 
Navier-Stokes equations, it is appropriate to examine the physical conditions that might be expected at 
the interface between liquid core and solid mantle. We shall therefore revert to the linearized version 
of the Navier-Stokes equations which was obtained in our previous work [1] and establish some funda¬ 
mental properties of its solutions. 
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For the scope of a first-order approximation, we can safely ignore the influence of the Earth’s 
rotation; however, we want to study those solutions of the Navier-Stokes equations that represent har¬ 
monic oscillations depending on the time factor exp ( i<rt ). We believe that a deeper understanding can 
be obtained of the conditions existing at the subject interface if we consider the Earth’s permanent 
deformations as the limit of spheroidal oscillations of the harmonic type when their frequency <r tends 
to zero. 

We use here the same notation adopted in the previous report; in particular, recall that the three 
sets of functions U ., K and R represent the coefficients of the spherical harmonic developments for the 
components u, v of the displacement vector, and for the potential due to deformation. 

Within the liquid core where p “ 0 and if we ignore the effects due to rotation, the equations 
governing the spheroidal oscillations of constant frequency cr can be obtained from Eqs. (39) of our 
previous report [ 1 ] in a greatly simplified form to be written as follows; 

tr 2 p 0 rV + p 0 R - yop a U + X* - 0; (43) 

<r 2 po U 4- pqR' + yopoX — po(yo^)’ + (XJf)' — 0; (44) 

r 2 R" + 2rR' — nOi + i )R - ^nGr 2 {p' Q U + p 0 X ). (45) 

For typographical convenience, we have omitted the subscript n referring to the order of the harmonic 
because our considerations henceforth will apply to any given value of n. 

In the preceding equations, primes denote derivatives with respect to the radial distance r, we 
have also set 

rX- rU'+ 2U- n(n + \)V (46) 

to represent the dilatation of the material. If we differentiate Eq. (43) with respect to r, subtract from 
it Eq. (44) and simplify the ensuing result by means of the original Eq. (43), we reach the following: 

U +~-\x-<T 2 (Y+rV’- U). (47) 

l Po I 

In the limit, as <r approaches zero while both U, V remain finite, we either reach the adiabatic condition 

yo + " 0, (48) 

Po 

also known as the Adams-Williamson condition, or we must have X — 0. This means that the dilata¬ 
tion must vanish. 

Pekeris and Accad have discussed the constitution of the liquid core by means of a nondimen- 
sional stratification function /3 (r) defined as: 

y° + - yo/90 - )- (49) 

(See Ref. 5). This in turn, is related to the Brunt-Vaisala frequency W J , according to 

fi - Kr N 1 . (50) 

Portf 

Excluding the existence of neutral stratification throughout the whole liquid core, i.e., /3 = 0 every¬ 
where, we realize that, as <r approaches zero X must also tend to zero. Soiving then Eqs. (43), (46) 
and (45) under these specific conditions, we find the particular solution 
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w - R*t y 0 ; V -r^TTT lr(u * r + 2U * ] ' 

n(n+l) 

r 2 (R*)" + 2 riR*) 1 - [n(n+l) + 4w Gpor 2 /y 0 ]R* - 0. 

We must choose that expression for R* that remains finite at the origin, r — 0; this will give rise to 
only one arbitrary parameter. Furthermore, if we allow that the variable V* be discontinuous at the 
core-mantle interface, we will have available only two parameters for the purpose of satisfying the three 
boundary conditions that are valid at the loading surface, r - a. 

One plausible physical condition that can justify the use of an extra free parameter is the existence 
of an infinitesimal boundary layer at the top of the liquid core, because within such layer the dilatation 
X can switch from zero to a finite value X*. This arbitrary jump in value can play the role of the addi¬ 
tional free parameter. 

We shall prove in a rigorous way that such a condition really exists at the core-mantle boundary; 
we shall accomplish this by considering the asymptotic expansion of the equations of motion with 
respect to the inverse frequency of oscillation. We shall limit our discussion to a simple physical model 
consisting of a uniform liquid core surrounded by a uniform solid mantle so that the analytical develop¬ 
ments can be simplified without detracting any essential characteristic feature from the phenomenon. 

Let us write the fundamental equations of motion for the case of a uniform liquid core, that is 
when po and A are constant and p. — 0. It is convenient to introduce the constant, 

. 4 _ 

A - y IT Op 0 , 

which is related to the gravitational acceleration according to 

yo(r) - Ar, 

and use it to define a new nondimensional constant, 

a — <r 2 /A, 

and a new variable Q — R/A, which has the dimension of a squared length. In this notation, Eqs. (43) 
and (45) become: 

~ X- rW- a V)- Q, (51) 

PqA 

r 2 Q"+ 2rQ’~ n(n + \)Q - 3r 2 X. (52) 

Equations (46) and (47) can be combined into a three-term relation: 

rX-a(V+ rV- U) - rU’ + 2U- n(n + 1)K (53) 

We must solve Eqs. (51), (52), and (53) for U , V, and Q. Let us begin with Eq. (53) and equate its 
second and third member; after some algebraic manipulations, we can write 

r(U - aV)'+ (2 + «)(l/-aF) - [n(n + 2) - a(a + 1)1 K 

It is easy to realize that r l+a is an integrating factor for the above equation; this will ultimately lead to 
its solution in the form 

U- aV- [n(n + 1) - a(a + l)]/(r), (54) 


where 


r 2 +°l{r)- fa r'+ a V(r)dr. 


( 55 ) 
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We next express the rX appearing in the right-hand side of Eq. (52) by means of the second member of 
Eq. (53), and eliminate the U through Eq. (54). The right-hand side of Eq. (52) can then be written as 

3ar{rV + (1 — a) V — [n(n + 1) — a(a + 1)]/). 

The preceding expression is instrumental in establishing the fact that 

3 arl(r) 

is a solution to Eq. (52). This can be verified by the use of the following relations, 

(r/)'= V - (1 +<*)/; 

r(rl)" - rV' - (1 + a) V + (1 + a)(2 + a)l\ (56) 

which are obtainable by differentiating Eq. (55) defining the integral Hr). 


The general solution of Eq. (52) can be written as 

0 - Br" + 3 arl(r). (57) 

where B is the only arbitrary constant appearing in our formulation because 0 must remain finite at 
r — 0. We shall use Eqs. (54) and (57) to eliminate the integral Hr) between them and express the 
dilatation X , as given by Eq. (51), in terms of 0 alone. This can be used to rewrite Eq. (52) as an 
equation containing only the unknown function 0. If we adopt the nondimensional independent vari¬ 
able 

s - r/a. 


the equation governing the variation of Q becomes 

s 2 - Q + 2s — n(n + 1)0 + —- [a(a + 4) — n(n + I)Js 2 0 — Ls" +2 . (58) 

ds l ds Ka 


Here L denotes the expression 

- [n(fl + 1) — «(a + \)]Ba n , 

Ka 

which depends on the arbitrary constant B. Equation (58) can be shown to have the particular solution 
Ds", where 

n fl(w + 1) — g(q + 0 an 
u ~ n(n + 1) — a(a +4) 

Its general solution is of the form 


0 - Ds" + r, 

where T is the general solution of the equation obtainable from Eq. (58) by deleting its right-hand side. 
The latter equation can be rearranged into 

■^r(sT)- IC(s)+p 2 ](sD-0, (59) 

ds 1 


where 


C(s) 


■ n(n -t 0 - (« + 4) £2^1, 

A 


and 


u 2 mm 


mini + 1 , . fsd^ 


n(n + 1) > 0. 


(60) 







PAOLO LANZANO 


When the frequency a-, and consequently the parameter a, decrease to zero and for s 0, the values 
of the function C(s) will become and remain negligible as compared with the constant v 2 . We 
therefore conclude that for 0 < s < s 0 , the asymptotic solution of Eq. (59) valid for decreasing values 
of <r is of the form 

r- - exp Ms - So)] + 0(a) 
s 

where F stands for an arbitrary constant and so — b/a, b being the radius of the liquid core. 

The asymptotic expansion of Q accordingly will be 

Q — Ds" + — exp Ms - s 0 )l + 0(a); (61) 

s 

and the asymptotic expansions of the other pertinent variables U, K and X can be evaluated therefrom. 


In fact, if we express the first of Eqs. (56) in terms of the nondimensional variable s, we get 

P(s) - s ~ + (2 + a)I(s). 
as 

The expression for /(s) is obtainable by equating Eq. (61) to Eq. (57) since they both represent the 
Q- variable: 

Q - Ds n + -j exp Ms - s 0 )] - Ba"s " + 3asa/(s). 

Proceeding along these lines and neglecting the positive powers of the small parameter a (or a), and 
since we are considering asymptotic expansions, we find, after elementary operations, that 

F(s) - —(as )" -1 + exp Ms - s 0 )]. 
n 3 aas 

Similarly, by using Eqs. (51) and (54), we shall reach the asymptotic representations of U and of the 
normal stress \X\ the results are: 

Uis) - Bias )" -1 + exp Ms — s 0 )l; 

3aa s i 

Xis) - - 5 - — exp Ms - s 0 )J. 

3 <r s 


If we introduce the nondimensional arbitrary parameters 

E\ - Ba n ~ 2 , F] - F/3aa 2 , 
we can write the previous results in dimensionless form as follows: 


- P(s) 

a 

- Uis) 
a 

Ris) 

ayoia) 

Xis) 


— E\ s" -1 + — i iF] exp Ms - s 0 )]; 
ft s 

£,s" -1 + —- - y- — F i exp Ms - s 0 )]; 
s l 

_ £<jL - [i - «-| £ + bfl exp Ms - s 0 )l; 

a 1 A a 2 ( n j s 


F i exp Ms - s 0 )]. 


(62) 


These equations prove that within the liquid core the oscillations induced by a long-period forcing func¬ 
tion vary as exp (1/cr) and not linearly with or 2 as one might have originally surmised by inspection of 
Eqs. (43) and (44). 
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The monomial terms that appear in the expressions for U and R are regular functions in the inter¬ 
val 0 < i < s 0 and will offset any discontinuity that might arise because of the presence of the 
exponentials. The F variable varies as v — Mo and will tend to increase; we know, however, from 
other sources that F will be discontinuous at the interface. 


On the other hand, the variable X related to both the dilatation and the normal stress in the liquid 
core depends on the exponential alone since <*v 2 is a finite constant; the X will therefore exhibit a 
discontinuous behavior typically induced by a boundary layer as described in the following. 

The behavior of the exponential term depends on the relative values of (s — s 0 ) and v Mo. 
For small values of o and for s much less than % s ~ s o >s negative, and the exponential of a large 
negative number provides a negligible contribution; we can say then that, from a practical point of 
view, the X variable vanishes at large distances from the top layer. 


On the other hand, consider the situation where s - s 0 is small enough to be comparable in mag¬ 
nitude with the assumed small value of o ; the exponential then provides a sizeable contribution. The 
smaller the value of o , the thinner the width of the topmost layer where X is nonvanishing. In the 
limit as a approaches zero, the stress distribution can be assumed to behave as a delta-function having 
a finite spike at $ — So- 


For s 


so 


b/a, Eq. (62) become 

- U(b) - E\ sS~' + - ( ~t 1 - F,; 
a stf 

X{b ) - — F, - n(n + 1) F,; 

So Xs 0 


R(b) 

ay 0 (a) 


E\ sfii 


(63) 


whereas Fean be taken to be discontinuous at the interface. We have two arbitrary parameters E h F, 
plus the value of the discontinuity suffered by Fin going from the liquid to the solid layer. These are 
then three quantities we can use in trying to satisfy the boundary conditions at r — a. 


This approach was followed in evaluating the load numbers of low order, e.g., of order n =■> 2, 3, 
4, which are rather well known from direct physical measurements. The agreement was good. We are 
therefore convinced that our analytical solution is an adequate representation of the physical conditions 
existing at the liquid core-solid mantle interface; and believe that the same procedure should be used 
for the calculation of higher order load numbers which are not known from ground data but which can, 
nevertheless, play a significant role in the final determination of the spheroidal deformations. 
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